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I. 



INTRODUCTION 



Surface ship and submarine survivability and cost have 
become the U. S. Navy's primary ship design concerns over the 
last several years. Unfortunately, these two criteria are 
often incompatible. In general, improved survivability leads 
to higher overall design and manufacturing costs for ships and 
submarines. In addition, congressionally mandated 

survivability testing has added substantial costs to the 
design and acceptance process for ships and submarines and 
their associated equipment. 

A key survivability issue is a vessels ability to 
withstand the shock and blast effects of an underwater 
explosion. This issue has been emphasized for both the 
surface and submarine communities by recent events. 
Operations in the Persian Gulf both during the Iran/Iraq war 
and during operations Desert Shield and Desert Storm required 
movement of ships in shallow, mine infested waters. There is 
little reason to believe that future conflicts will be any 
different . 

Currently recommended survivability goals state that it is 
desirable that ships should be able to stay afloat and 
continue to fight after a hit by enemy ordinance similar to 
the mines encountered by the USS Roberts, USS Tripoly and USS 
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Princeton. The near loss of the USS Roberts emphasized a need 
to improve the structural strength of warships. While the 
Princeton incident clearly demonstrates the positive impact of 
the Navy shock hardening program, it also shows that the Navy 
has not yet completely achieved the stated survivability goal. 

Similarly, especially with a reduced deep water threat 
resulting from the dissolution of the Soviet Union and the 
expected division of the Soviet Navy, U. S. Navy submarines 
deployed in future conflicts are more likely to be required to 
operate against small conventionally powered submarines in 
shallow coastal waters with the same mine threat already 
encountered by surface units. 

Therefore significant interest continues in determining 
methods that can reduce costs in both the design and 
manufacturing phases of ship procurement while still improving 
the shock and structural hardness of both surface vessels and 
submarines. One primary means of achieving cost savings 
without loss of desired attributes is the process of 
optimization of ship designs. This process reduces cost by 
removing unneeded redundancy. However, to perform an 
optimization, a designer must have a clear understanding of 
the modes of failure and their causes. Unfortunately, the 
interaction between an underwater shock wave and a ship hull 
as well as the effect on equipment is a highly complex event 
and remains poorly understood. In the past, data for studying 
failure modes and the underlying physics of the events leading 
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to failure has been obtained through costly and time consuming 
process of underwater explosion testing. This process, 
although useful, is extremely limited. Models tested to 
failure can only be used once and each model can only contain 
a limited number of data collection instruments. In addition, 
environmental safety concerns are increasing costs and 
imposing limitations on the Navy's ability to carry out 
required testing. 

The rapid improvement of digital computing technology over 
the last decade is a possible solution to the problem. While 
computer technology cannot completely replace live testing, 
properly constructed computer models provide an opportunity to 
obtain an unlimited amount of data without the limitations 
associated with live testing. In addition, if computer model 
predictions can be shown to be reliable, some cost reduction 
may be possible through the reduction of live testing required 
to meet congressionally mandated survivability testing. 
However, currently existing computer codes used to construct 
and process computer models in this field, although promising 
are not yet proven . 

A research program is underway at the Naval Postgraduate 
School to study numerical modeling of ship structures 
subjected to both near and far field underwater explosions. 
This program is expected to improve the understanding of 
factors affecting the reliability of numerical models and 
provide insight into the dynamic response of surface ship and 
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submarine hulls and the physics that lead to failure when a 
hull is subjected to an underwater shock wave. The current 
study centers around simple cylinders constructed of a 
homogenous material . Previous work in the program is 
documented in reference 1. Future studies will include more 
complex materials and structures as experience increases and 
the reliability of the numerical models is proven. 

This thesis deals specifically with the non-linear 
response of a cylindrical shell subjected to both near field 
and far field side-on attacks. The first object of the 
research was to prove that the numerical computer code 
developed for the research provided correct results. The 
verification included comparison of numerical results to 
analytical solutions for two simple models. The results show 
that the numerical solution closely matched the analytic 
solution in both cases. 

The near field study compares numerical predictions with 
expected physical results in the shell and the surrounding 
water. It also discusses physical findings related to the 
near field attack. Specifically, the most significant damage 
was found to be local in nature. Damage included severe local 
buckling of the shell at the point nearest the charge. 
Further, stiffeners in stiffened models failed by local 
tripping . 

The far field study compares numerical predictions with 
experimental results obtained from a underwater explosion test 
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of an aluminum cylinder subjected to a far field side on 
attack. In addition, analyses were performed to determine the 
sensitivity of the results to: 

- mesh refinement 

- boundary effects 

- rotation from expected configuration 

- number of quadrature points used for integration 

- time integration increment 

- type of shell element formulation 

used for the analysis. Results show that numerical computer 
codes can generally match experimental results if end effects 
and rotation are correctly modelled. Inconsistencies between 
experimental and numerical results is most likely caused by 
uncontrolled factors associated with the underwater explosion 
test and overall element averaging rather than a failure of 
the numerical method to provide correct results. 

Physical results included observation of the primary 
response and damage modes. Primary response modes included 
accordion motion, whipping motion and breathing motion in the 
direction perpendicular to the charge. Damage was found to be 
more global in nature than that found for the near field case. 
Primary damage areas were located near the end plates on the 
side of the cylinder nearest the charge and at the center of 
the cylinder on the side most remote from the charge. 

Recommendations will be provided to improve control of 
future tests to improve final results. Finally, preparations 
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for future testing will be described and recommendations for 
additional study are provided. 
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II. NUMERICAL COMPUTER CODE 



A. GENERAL DESCRIPTION OF THE METHOD 

The near field portion of this study was performed using 
the finite element method. The finite element method uses a 
discretized representation of the cylindrical model and the 
surrounding water along with appropriate material models and 
equations of state to obtain a solution for the response of 
the shell to a near field explosion. 

The far field study was performed using the finite element 
method to obtain the response of the shell. However, the 
surrounding water was modelled using the boundary element 
method. The boundary element method uses fundamental physical 
principles to reduce the surrounding media and the associated 
forces to discrete forces and masses which are applied to the 
nodes of a two dimensional mesh. The two dimensional mesh, in 
turn is superimposed over the surface of the shell. Use of 
the boundary element method significantly reduces the number 
of elements required for the numerical model, and subsequently 
reduces the computational effort required to solve the 
problem. For example, models used for the near field problem 
consisted of more than 13000 three dimensional elements. In 
contrast, the most refined model used for the far field 
problem consisted of only 984 elements, or approximately 7.5% 
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of the total required for the near field problem. In addition 
to the reduced computational effort, the reduced number of 
elements allowed the storage of significantly more data per 
element, providing additional flexibility and improvement in 
the study. The finite element method (FEM) computer code 
used for the analysis was VEC/DYNA3D and the boundary element 
method (BEM) computer code was USA (Underwater Shock Analyzer) 
code . 

To solve the far field problem, a link was developed to 
allow the boundary element computer code and the finite 
element computer code to operate interactively. This linkage 
was developed in 1991 at the request of the Naval Postgraduate 
School under funding provided by the Defense Nuclear Agency 
( DNA) . A more detailed description of the two computer codes 
is provided in the following sections. 

B. VEC/DYNA3D FINITE ELEMENT METHOD CODE 

VEC/DYNA3D [Ref. 2] is an explicit finite element code. 
It has been used successfully for various types of dynamic 
nonlinear engineering problems since its conception in 1976. 
VEC/DYNA3D was selected for this study for several reasons. 
First, as stated above, VEC/DYNA3D is an explicit code. This 
attribute has two distinct advantages and two disadvantages. 
The advantages are its relatively high speed and its ability 
to be implemented on a relatively small stand alone 
engineering workstation. Initial work for this study is being 
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performed on IBM RISC 6000 workstations. Once the USA/DYNA3D 
interface is proven to be reliable and accurate and experience 
has been gained in the use of the software, work will begin on 
more complex models using main frame type computers. 
Therefore it was important to obtain a code that was able to 
work significant problems on a small workstation and yet be 
compatible with the main frames expected to be used in the 
future. VEC/DYNA3D is compatible with a full range of 
engineering workstations and has been implemented on the Los 
Alamos CRAY computer. Problems including up to 20000 solid 
elements have been run on Naval Postgraduate School 
workstations with 16 megabytes of random access memory. 

The first disadvantage associated with the explicit 
numerical code is that the code is not inherently stable. 
This means that any problems dealing with time integration, 
including the underwater shock problems included in this study 
must be treated with care. Integration time steps must be 
matched closely with the size of the elements in the problem. 
This is performed automatically by VEC/DYNA3D in the stand 
alone mode. However, when coupled with the USA code, this 
automation is no longer functional. Incorrect selection of 
integration time steps can lead to significant oscillations 
and/or inaccuracies in the final solution. 

The second problem associated with the explicit codes is 
the mesh reflection effect. In explicit codes, non-uniform 
meshes result in inaccurate solutions due to mesh reflection. 
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Two factors appear to be important in ensuring that correct 
solution was obtained. The first is mesh size and the second 
is total mass of neighboring elements. Sensitivity analyses 
indicate that error in the final solution is relatively small 
if neighboring elements are kept within ten percent of each 
other in size. This rule of thumb was used whenever possible 
in the performance of this study. However, this factor leads 
to some inefficiency in obtaining solutions since refinement 
must be performed over a larger area of the mesh to obtain a 
mesh independent solution than might normally be required in 
an implicit code. The additional area means more total 
elements and a subsequent increase in computation time to 
obtain the problem solution. Careful selection of the time 
integration increment can minimize this problem and allows use 
of larger variations at the expense of longer problem run 
times. These disadvantages can be overcome through careful 
planning. In general, they did not significantly overshadow 
the benefits associated with using an explicit code. 

The second reason for selecting VEC/DYNA3D was its high 
degree of flexibility. It has a wide range of available 
material models and equations of state including the ability 
to model strain rate sensitivity, explosive materials and 
acoustic media. This flexibility is enhanced by the large 
degree of interactivity available with VEC/DYNA3D when used 
with the INGRID pre-processor [Ref. 3] and TAURUS post- 
processor [Ref. 4] . Changes can be entered with relative ease 
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using the pre-processor and almost any desired quantity can be 
obtained through knowledgeable use of the post-processor once 
the calculations are complete. 

The final reason for selecting VEC/DYNA3D is that it 
is available in the public domain. Therefore the only cost 
associated with its use in this project was the cost 
associated with developing the link between VEC/DYNA3D and 
USA. 

C. USA BOUNDARY ELEMENT METHOD CODE 

The Underwater Shock Analyzer (USA) code [Ref. 5] is a 
boundary element code based on the Doubly Asymptotic 
Approximation (DAA) theory developed by T. L. Geers in 1971 
[Ref. 6] . Through the use of the DAA theory and the boundary 
element formulation, USA computes the acoustic pressure 
loading and added mass matrices which represent the fluid 
surrounding the submerged shell. The acoustic pressure 
loading and added mass are applied at selected wetted nodes. 
This formulation has the benefit of significantly decreasing 
the number of elements required to model the submerged system 
since external water elements need not be included in the 
calculations. As stated previously, the reduced number of 
elements requires substantially less time and storage space to 
obtain a solution. 

A detailed description of the DAA theory used in the USA 
code is provided by reference 7. In general, the governing 
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differential equation for the structure and the acoustic wave 
equation along with geometric compatibility between the fluid 
and the structure and appropriate initial conditions and 
boundary conditions are used to solve for the general 
structural response. 

The governing differential equation for the structure is 
given below. 



[M s ] {x} + [C s ] {X} + [K s ] {X} = {f} (1) 

where [M s ] , [C s ] and [K s ] are the structural mass, structural 
damping and structural stiffness matrices. Further, 

{x} , {x} and {x} are the vectors of nodal acceleration, 
velocity and displacement for the structure. 

The {f} vector is the vector of nodal excitation forces 
generated by the acoustic wave. The nodal excitation force is 
a function of the incident and scattered pressures of the 
impinging acoustic wave and any concentrated loads applied to 
the structure. The equation specifying this relationship is 
shown below. 



{f} = -IG) [A f ] {P, + P s } + {f d } (2) 

where [G] is the fluid/structure transformation matrix, [A f ] 
is the diagonal area matrix associated with the fluid mesh, 
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{f d } is the vector of concentrated forces and {P i + P s } is the 
vector of incident and scatter pressures. 

In the above equation, all quantities are known with the 
exception of {P s } . The scattered pressure can be represented 
in the form of a first order ordinary differential equation in 
terms of the fluid mass, element areas, density of the water, 
speed of sound in water and the vector velocity of the 
scattered acoustic wave. 

[M f ] {P s } + p c[A f ] {P s } = p c[M f ] { u s } (3) 

In the above equation p, c, [M f ] and {u s } are the water 
density, speed of sound in water, the fluid mass matrix and 
the vector of acoustic velocities respectively. The remaining 
quantities were previously defined. The fluid mass matrix is 
a diagonal matrix of virtual masses calculated by the boundary 
element method. This virtual mass is added to the structural 
mass and becomes important at low structural frequencies where 
the inertia of the water surrounding the structure can 
significantly affect the response of the structure. 

At this point, the doubly asymptotic approximation is 
applied. The first asymptotic approximation, or high 
frequency approximation applies at early times very near the 
surface of the structure. It assumes that |P S | > |P S | . This is 
a logical assumption and can be easily visualized in the one 



13 



dimensional case since the scatter wave of a one dimensional 



wave will completely reverse its direction according to 
Snell's law when it strikes a nearly rigid structure. As the 
change in direction is virtually instantaneous, the time 
derivative of {P s } , or in other words {P s } , will be nearly 
infinite. Therefore, after integration, the differential 
equation which describes the scattered wave pressure is 
reduced to: 



{P s } = p c{u s } (4) 

allowing the direct solution of the scattered wave pressure. 

The late time asymptotic approximation, or low frequency 
approximation assumes that |p s | < |p | . Again, this is a 
logical assumption since the velocity of the scattered wave 
becomes constant at the speed of sound in water as the 
scattered wave travels a significant distance from the 
structure surface. As a result, the differential equation 
describing {P s } can be modified as shown below. 

[A f ] { P s } = [M f ] {u s } (5) 



The above solution is known as the first order Doubly 
Asymptotic Approximation (DAAl) and is exact at early and late 
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times for a planar surface. However, it does not take into 
account the effects of curvature. 

Since most surfaces have some curvature, a second DAA 
theory was developed called DAA2 or the second order Doubly 
Asymptotic Approximation. The DAA2 theory is described in 
reference 8. 

It must be noted that this code has limitations which 
result directly from the fundamental assumptions associated 
with the DAA theory [Ref. 5] . First, DAA is not theoretically 
appropriate for concave or multiple structures or near surface 
problems involving convex bodies. However, studies show that 
only results in highly shadowed, closely spaced areas or 
regions of strong concavity are affected. Secondly, DAA 
requires that the source of the incident wave be sufficiently 
removed from the structure since it can only account for 
acoustic waves and not hydrodynamic flow. Finally, the DAA 
theory is based on an early time (high frequency) 
approximation coupled with a late time (low frequency) 
approximation. Therefore, although the DAA solution will be 
very good at early times when the high frequency approximation 
is dominant and at late times when the low frequency 
approximation is dominant, it can vary significantly from the 
analytic or exact solution during intermediate times when 
neither the high or low frequency solution is dominant. 
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D. USA/DYNA3D CODE VERIFICATION 



Since the USA/DYNA3D interface was new and had not been 
tested, some effort was expended on performing a verification 
of the performance of the combined computer code. To perform 
the verification, two cases with known analytic results were 
modelled using the USA/DYNA3D code. The first case was a 
quarter sphere and the second was an infinite cylinder. 
Results were satisfactory for both cases and the code 
interface is believed to be performing correctly. 

1. DETAILED DESCRIPTION OF THE SPHERICAL MODEL 

Figure II. 1 shows the quarter symmetry model used for 
comparison to analytical results. The model contains 150 





Figure II. 1. Elastic sphere verification model. 
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elements. Figure II. 2 shows the verification model test 
geometry. The thickness to diameter ratio of the shell is 1 
to 50 and the shell is constructed of steel. The excitation 
is provided by a very small step pressure wave. As a result. 




Figure II. 2. Elastic sphere verification model geometry. 

the shell response is considered to be completely elastic. 
The case was run using the elastic material model of DYNA3D 
and, since results are being compared to the analytic results 
found in reference 7, the same material and water properties 
as those found in reference 7 were used. As stated in 
reference 7, the exact solution is obtained from separation 
of variables as shown in reference 10. The material and water 



17 



properties used are listed below: 



E=206 . 84 GPa 
V =0.33 

p =7784.5 kg/m 



Steel Properties 
Young's Modulus 
Poisson's ratio 
Mass density- 
water Properties 

Sound speed c=1461.2 m/s 

Density p=999.6 kg/m 

The numerical results using the USA/DYNA3D combination 
for the above test case compare favorably with the exact 
results. The normalized results are shown in Figure II. 3. it 
can be seen that the numerical results lag the exact results, 
but the difference is negligible. 

2 . DETAILED DESCRIPTION OF THE INFINITE CYLINDER MODEL 



The infinite cylinder model was run using the same 
material and water properties shown above. Figure II. 4 shows 
the geometry used for the analyses and, as shown in Figure 
II. 5, a single ring of elements was used to model the infinite 
cylinder by enforcing symmetry boundary conditions on each end 
of the model. In addition, since this is a two dimensional 
problem, the TWODIM option in USA was used to generate the 
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Figure II. 4. Infinite cylinder verification model. 
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Figure II. 5. Infinite cylinder verification geometry. 



added mass and DAA matrices. Further, the value of the T| 
variable was set to 0 . 0 . T) is a factor that accounts for 
curvature in the structure. If T| = 0.0, there is no 
correction and the DAA solution provided is equivalent to the 
DAA1 solution which is exact for a planar surface. T| - 1.0 is 
suitable for a spherical surface. Structures that are neither 
planar or spherical require a value between 0.0 and 1.0. In 
this case, the use of T] = 0.0 allowed computation of a DAAl 
solution for comparison to a known analytic DAAl solution for 
a cylinder. The first model attempted had a longitudinal 
length of 0.1 inches. However, it was discovered that this 
resulted in a oscillatory solution as shown in the first graph 
in Figure II. 6. A similar oscillation occurred on the reverse 
side of the cylinder as shown in the first graph of Figure 
II. 7. After a check of the input data to ensure that the 
problem was not caused by numerical instability, it was 
hypothesized that the oscillation was caused by residual three 
dimensional effects caused by the finite width of the model. 
As a test of this hypothesis, two additional models were run 
with widths of 0.01 and 0.001 inches. As shown in Figure II. 6 
and II. 7, reduction in width progressively reduced the 
oscillations on both the front and back of the cylinder. At 
0.001 inches, oscillations are absent. 
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NORMALIZED VELOCITY (pcu/P,) NORMALIZED VELOCITY (pcu/P,) NORMALIZED VELOCITY (pcu/P,) 



L/a = 0.1 . n =0.0 




NORMALIZED TIME (ct/a) 
L/a = 0.01, n = 0.0 




NORMALIZED TIME (ct/a) 
L/a = 0.001, 7i = 0.0 




NORMALIZED TIME (ct/a) 



Figure II. 6. Near element oscillation reduction for 
infinite cylinder verification model. 
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NORMALIZED VELOCITY (pcu/P,) NORMALIZED VELOCITY (pcu/P,) NORMALIZED VELOCITY (pcu/P,) 



L/a = 0.1 , t| =0.0 




NORMALIZED TIME (d/a) 
L/a =0.01 , n = 0.0 




NORMALIZED TIME (a/a) 
L/a = 0 001, ti = 0.0 




NORMALIZED TIME (a/a) 



Figure II. 7. Remote element oscillation reduction for 
infinite cylinder verification model. 
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The final results from the 0.001 inch model with an 



T| variable value of 0.0 were compared to the analytical exact 
and analytical DAA1 solutions with favorable results as shown 
in Figure II. 8. It can be seen that the results on both the 
front and back sides of the cylinder lie very close to the 
analytic DAAl solution. 

A further investigation was conducted to determine 
what value of the T| variable would result in the numerical 
solution closest to the analytic modal solution. Values of 
0.0, 0.25, 0.5, 0.75 and 1.0 were tested. Review of the 
results show that the value of the T) variable that provides 
the results nearest the analytic modal solution varies 
depending on time and position on the cylinder. 

For the front of the cylinder, an r| variable value of 
0.0, as shown in the first plot Figure II. 8, provides results 
fairly close to both the analytic modal and analytic DAAl 
solution for values of normalized time less than four and one 
half. However, as would be expected with an T) value of 0.0, 
the numerical solution approaches the analytic DAAl solution 
at late times. Figures II. 9 and II. 10 show the results for T) 
values of 0.25 and 0.5. Although they do not match either the 

t 

analytic modal or analytic DAAl solutions at early times, they 
do come closer to matching the analytic modal solution at late 
times with T) = 0 . 5 providing the best results. 
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NORMALIZED VELOCITY (pcu/P,) 



NEAR ELEMENT 




NORMALIZED TIME (ct/a) 



REMOTE ELEMENT 




Figure II. 8. Infinite cylinder results with r| = 0.00. 
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NORMALIZED VELOCITY (pcu/P,) 



NEAR ELEMENT 




REMOTE ELEMENT 




NORMALIZED TIME (ct/a) 



Figure II. 9. Infinite cylinder results with T| = 0.25. 
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NORMALIZED VELOCITY (pcu/P,) 



NEAR ELEMENT 




NORMALIZED TIME (ct/a) 



REMOTE ELEMENT 




NORMALIZED TIME (ct/a) 



Figure II. 10. Infinite cylinder results with T| = 0.50. 
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On the reverse side of the cylinder, a value of 0.0 
provides a result very near the analytical DAAl solution. 
Again, this is the expected result with a T| variable value of 
0.0. However, as can be observed in the second plot of Figure 
II. 8, the numerical solution varies substantially from the 
analytic modal solution between the normalized time values of 

I. 0 and 4.0. Values of 0.5, 0.75 and 1.0, as shown in Figures 

II. 10, II. 11 and 11.12, provide results near the analytical 
modal solution with 0.75 coming the closest. 

Assuming that interest lies in late time results over 
the entire cylinder, the results show that the best overall 
value of r| for an infinite cylinder lies between 0.5 and 0.75. 
More compact bodies will have best results with higher T| 
values. Based on the above results, an r| value of 0.5 was 
selected for use for most models used for analyzing the far 
field problem described in Chapter IV. 
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NEAR ELEMENT 




NORMALIZED TIME (ct/a) 
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Figure II. 11. Infinite cylinder results with T| = 0.75. 
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NORMALIZED VELOCITY (pcu/P,) 



NEAR ELEMENT 




NORMALIZED TIME (ct/a) 



REMOTE ELEMENT 




Figure 11.12. Infinite cylinder results with T) = 1.00. 
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III. NEAR FIELD SIDE-ON UNDERWATER EXPLOSION ANALYSIS 



A. GENERAL DESCRIPTION AND OBJECTIVES. 

The first case studied using the previously described 
numerical tools was a simple cylinder subjected to a near 
field side-on explosion. The study included both stiffened 
and unstiffened cylinders. The main objective of this study 
was to further the understanding of phenomena associated with 
a near field explosion. Of specific interest were shock wave 
propagation, water/shell interaction and nonlinear response of 
the shell and stiffeners. A clear understanding of these 
phenomena is required to determine the damage mechanisms that 
lead to shell collapse. Understanding these factors would 
ultimately allow the optimization of underwater structures to 
withstand this type of attack. Since there was no 
experimental data for comparison, the accuracy of the model is 
based on observing certain key expected performance 
parameters. Future studies should include experimental to 
numerical comparisons to validate the model and its results. 
The following sections provide a detailed explanation of the 
models used and provide results obtained from the study. Much 
of this information has already been published in reference 11. 
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B. DESCRIPTION OF MODEL USED FOR THE STUDY. 

The model used for this study consisted of an explosive 
charge, a cylindrical shell and the surrounding water. Since 
the problem setup included a symmetric structure with 
symmetric loading, quarter symmetry analysis was used to 
minimize computational effort while providing the greatest 
opportunity for refinement of the discretization mesh. 
Different modeling techniques were used for each section of 
the model. Figure III.l shows the combined quarter symmetry 
model used for the analysis. 




Figure III.l. Quarter symmetry model for near field side- 

on attack analysis. 
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1. CYLINDER MODEL 



The quarter symmetry discretization of the cylindrical 
shell is shown in Figure III. 2. The discrete model represents 
a physical cylinder with a 42 inch length and 12 inch 
diameter. The axial surface is one quarter inch thick and 
each end plate is has a one inch thickness. For the stiffened 
model, two one eighth thick by one inch high circumferential 
stiffeners were spaced equidistantly within the shell. 




Figure III. 2. 



Shell quarter symmetry discretization for 
near field side-on attack analysis. 
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The material selected for the analysis had a yield 
strength of 43 KSI and a Young's modulus of 10800 KSI . An 
elastic-plastic constitutive model with no strain or strain- 
rate hardening was selected. The constitutive model and 
material properties are consistent with 6061-T6 aluminum. The 
end plates are considered to be rigidly attached to the shell. 

The Belytschko/Lin/Tsay element formulation was used 
for the model . This formulation was selected for its 
computational efficiency and high degree stability in the 
presence of large deformations. The theory associated with 
this formulation is discussed in reference 12. 



2 . EXPLOSIVE CHARGE 

The shock wave for the study was provided by an 
explosive charge placed with its center located one foot below 
the geometric centroid of the cylinder. This resulted in a 
six inch standoff distance between the center of the charge 
and the surface of the model . 

The Jones /Wilkins/Lee (JWL) equation of state was used 
to model a one inch radius pentolite charge for the 
unstiffened study. A slightly larger charge with a radius of 
1.2 inches was used for the stiffened model. The JWL equation 
of state as stated below: 



p-Ail - ad - -^L) 



R,V 



r 2 v 
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determines the pressure generated by the charge. It has an 
empirical basis and is commonly used to model the detonation 
products of high explosives. The coefficients A, B, C, 
R 1 , R 2 and co are experimentally derived and tabulated for each 
type of explosive. Additional information concerning the JWL 
equation of state can be found in reference 13 . 

Since the explosive charge has a physical volume, it 
is consumed at a rate governed by its burn rate and 
distribution about the detonation point. The combination of 
the equation of state with the finite burn time generates a 
shock wave in the surrounding media that is a function of 
explosive material weight, explosive material distribution, 
time, and distance from the charge. The quarter symmetry 
discretization of the explosive charge is shown in Figure 
III. 3. The quarter charge is composed of 1296 solid elements. 
The extremely fine discretization was necessary to ensure 
spherical propagation of the shock wave while still 
maintaining a coherent mesh as the explosive material expanded 
in the surrounding water. 

3. SURROUNDING WATER 

As stated earlier, the boundary element method is not 
considered to be suitable for modelling for the near field 
case because of the fundamental theory that forms the basis of 
the doubly asymptotic approximation. As a result, the 
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Figure III. 3. Explosive charge quarter symmetry 

discretization for near field analysis. 

surrounding water medium was discretized using the finite 
element method. 

The water media discretization can be seen in Figure 
III.l. Even though quarter symmetry was used, the mesh still 
consists of 11640 solid elements. This proved to be a 
significant limitation in the use of the finite element method 
for the study of the underwater explosion phenomena. This 
limitation occurs because any increase in distance between the 
cylinder and the charge results in a significant penalty in 
the form of additional computational time and storage 
requirements as the size of the water medium was expanded to 
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accommodate the charge and cylinder. This problem can be 
overcome by the use of brute force - i.e. the use of larger 
and faster computers. However, this results in significant 
increases in cost in the performance of required studies. 
Fortunately, the boundary element method discussed in chapter 
II can be used at distances only slightly greater than the one 
foot distance used for this study. 

The water was modelled as an infinite medium. Since 
an infinite medium cannot be discretized, the infinite case 
was modelled using a finite volume bounded by non-reflecting 
boundaries. The non-reflecting boundaries cause a travelling 
wave to damp out so that there is no reflection from the 
boundary surface. The non-reflecting boundary was used on 
four sides of the water volume. The other two sides were 
symmetric boundaries since quarter symmetry was used in the 
problem . 

The Gruneisen equation of state was used to model the 
surrounding water. Since water cannot withstand a tensile 
force, a slightly positive pressure cut off was used in 
addition to the Gruneisen state equation to determine when 
cavitation occurred in the fluid. The pressure cutoff 
prevents the occurrence of physically impossible negative 
pressures in the water media. 
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C. RESULTS OF THE NEAR FIELD STUDY 

1. FREE FIELD TRANSMISSION AND SHELL WATER INTERACTION 

Figure III. 4 shows the water mesh located between the 
charge and the nearest point of the cylinder to the shell. 




Figure III. 4. Location of specified elements for near 

field analysis. 

The pressure histories for elements designated A through F in 
Figure III. 4 are plotted in Figure III. 5. Unless otherwise 
stated, time scales are in microseconds and pressure and 
stress scales are in grams per c/n-nsec 2 • Of the water 
elements shown in Figure III. 4, element A is nearest the charge 
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Figure III. 5. Pressure-time history for specified water 

elements for near field analysis. 
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and element "F" lies next to the structure. Several 
observations can be made on the basis of the Figure III. 5 
plot. First, the free field peak pressure (elements "A" 
through "E") drops as distance from the charge increases. The 
rate of decrease is slightly less than the 1/R reduction 
expected for a spherical charge. This deviation can be 
attributed to the nearness of the charge. Element "A" is 
between two and three charge radii from the charge while 
element "E" lies between five and six charge radii from the 
charge. Normally, the 1/R performance is expected to hold 
outside ten charge radii from the charge. Nearer distances 
are expected to decrease at less than the 1/R rate. 
Therefore, the pressure decrease with respect to distance from 
the charge is considered to be physically correct for this 
situation . 

The second observation is that the pressure suddenly 
increases in element "F" . This is an expected phenomena. By 
using Snell's law, it can be shown that pressure will double 
at the water/structure interface when an acoustic wave 
impinges on a rigid structure. In this case, the pressure did 
not completely double. The reduction in the peak can be 
attributed to two factors. First, the pressure provided for 
the element is an average of the entire element. In theory, 
the doubling only occurs at the water/structure interface. 
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Since the pressure in the other areas of the element are lower 
than at the actual interface, the average for the element will 
be lower than the interface pressure. Second, the structure 
is not completely rigid as assumed by the theory that predicts 
the doubling. As a result, the reflected pressure will be 
somewhat less than double the pressure of the impinging 
acoustic shock wave. 

It can be further observed that the pressure in 
element "F" drops rapidly to zero after peaking. This occurs 
because the impact of the shock wave on the shell surface 
causes the shell to deform inward at a speed higher than the 
adjacent water. Since water is highly incompressible, this 
velocity differential causes a rapid drop in pressure which in 
turn leads to a significant area of local cavitation. Figures 
III. 6 and III. 7 show the sudden pressure rise caused by the 
shock wave reflection and subsequent cavitation in the form of 
fringe plots. In the fringe plots, the darker areas represent 
areas of higher pressure and similarly, the lighter areas 
represent areas of low pressure. 

Finally, it can also be observed in Figure III. 5, that 
the reflected shock wave travels back toward the originating 
explosive charge. As stated previously, the reflected peak is 
visible in "F" as the rapid increase in pressure. In element 
"E", it is shown by the flattening and then continued drop of 
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Figure III. 6. 



Pressure fringe plot as the pressure wave 
strikes the nearest shell surface. 
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Figure III. 7. Pressure fringe plot as cavitation occurs 

in water near the shell surface. 
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the pressure plot starting at 75 microseconds. The cavitation 
next to the shell prevents the pressure in element "E" from 
developing into a peak. The reflected peak is clearly visible 
in the element "D" plot at 90 microseconds and can be seen to 
be developing first in "B" and then in "A” at 100 
microseconds. There are several key facts to note here. 
First, the height of the reflected peaks decreases as distance 
increases from the point of reflection. Second, the peaks 
show up in elements near the shell before they show up in 
elements further from the shell. In other words, as would be 
expected, elements are affected in the reverse order that they 
are affected by the incident shock wave. Element "C" appears 
to be an anomaly to this trend since no reflected pressure 
peak develops for "C" between the peaks for "D" and "B" . 
However, careful observation shows that the leveling off of 
the "C" element pressure plot is caused by the arrival of the 
rising reflected wave at the same time as the second peak of 
the incident pressure wave is decreasing. The result is a 
cancellation effect, and the reflection peak does not develop. 

Several other observations were made with regard to 
the free field pressure propagation and water/shell 
interaction. First, the free field pressure dissipates very 
quickly. By the time the shock wave has expanded to a point 
halfway up the side of the shell, the peak pressure has 
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already dropped to a value that is 20% of the free field 
pressure at element "E" . By the time it reaches the most 
remote point on the shell midway between the two ends, the 
peak pressure has further dropped to seven percent of the peak 
value in element "E" . These values are consistent with the 
1/R thumbrule for peak pressure reduction for a spherical 
charge. The value at the remote side is a little low, but 
this can be attributed to shadowing by the cylinder. 

2. SHELL RESPONSE 

a. UNSTIFFENED MODEL 

Figures III. 8 and III. 9 are fringe plots showing 
the propagation of stress in the shell shortly after initial 
impact of the shock wave with the shell. It can be observed 
that the stress propagation is initially a quarter circular 
shape. This is the result of quarter symmetry with elastic 
deformation. However, as plastic deformation begins in the 
shell, the stress pattern develops into the mushroom shape 
shown in Figure III. 9. This mushroom shape is caused by a 
combination of two effects. First, the shock wave is causing 
the cylinder to bend globally at the center like a beam. This 
tends to put the whole thickness of the surface area nearest 
to the charge into compression. This effect decreases as the 
vertical distance between the location on the shell and the 
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Figure III. 8. Propagation of Von Mises stress at 90 

microseconds (30 microseconds after 
pressure wave impact) . 
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Figure III. 9. Mushroom shaped propagation of Von Mises 

stress at 150 microseconds. 



46 




neutral plane of the cylinder decreases. It also decreases as 
the location gets closer to the ends since the maximum moment 
will occur at the centerline between the two end plates. At 
same time, the shock wave is deforming the surface in a 
circular pattern. This tends to relieve the compressive 
stress along the bottom of the cylinder in the vicinity of the 
shock wave impact. The same deformation tends to reinforce 
the compressive strain along a circumferential line running 
through the near point perpendicular to the axis. This in 
turn causes the mushroom shaped pattern observed in Figure 
III. 9. As the pressure wave continues along the surface of 
the cylinder nearest to the charge and circumferentially 
around the cylinder surface, the mushroom grows. As the shock 
wave passes the end of the cylinder, the pattern resumes its 
former quarter circular shape. 

Although the stress pattern changes shape, the 
effective plastic strain changes only marginally from a 
quarter circular shape to a quarter elliptical shape with the 
major axis in the same direction as the axis of the cylinder. 
This is consistent with the above description since the strain 

pattern covers a relatively small area of the lower part of 

« 

the cylinder and tends to spread after the mushroom stem has 
already passed the deforming location. 
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Since the speed of sound in aluminum is higher 
than in water, the stress wave propagated through the cylinder 
faster than through the water. As a result, it took only 120 
microseconds for the stress wave in the aluminum to arrive at 
the point on the mid shell circumference most remote from the 
charge while it took 200 microseconds for the wave to 
propagate to the same point through the water medium. Based 
on a linear travel distance of one foot and an acoustic 
velocity of 5000 ft/sec in water, the 200 microseconds appears 
correct. A similar check of expected travel time through 
aluminum based on an acoustic velocity of 15000 ft/sec with a 
circumferential travel distance of 1.5707 feet indicates that 
the 120 microseconds is slightly longer than the expected 105 
microseconds. However, plot states were taken only every 20 
microseconds. Therefore, if arrival occurred just after the 
previous plot state was recorded, it is possible for timing to 
be off by as much as 20 microseconds. The 15 microsecond 
difference noted above falls well under the 20 microsecond 
plot state difference. 

Mid plane effective (Von Mises) stress and 
effective plastic strain were plotted for selected locations. 
Effective plastic strain is defined by the relation: 

E p = ^ 2 “ ^ ^ e Pl ~ E P2 ^ 2 + ^ C P2 ~ e p3 ) 2 + ( e p3 ~ e pl ) 2 ] 2 
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with e pl , e p2 and e p3 representing the true plastic strain 
components [Ref. 9]. Figures III. 10 and III. 11 show the 
results for the shell elements nearest the charge. Based on 
comparison of the two plots, it is apparent that it took 
approximately 30 to 40 microseconds for the stress to build to 
the yield stress level. Once yield stress was attained, 
stress did not change since the perfectly plastic material 
model with no strain or strain rate hardening was used for the 
analysis. It is also clear from comparison of the two plots 
that plastic deformation continued as long as the material 
remained at the yield stress level. Plastic deformation 
discontinued and remained constant when the stress dropped 
below the yield level . Since the strain plotted was effective 
plastic strain, no strain reduction or unloading occurred as 
the material stress was reduced. It is also apparent from the 
effective strain plot that the initial loading period is the 
most important. Even though the stress remained at the yield 
stress for most of the first 1300 microseconds following 
initial impact of the shock wave, over 80 percent of the 
strain occurred in the first 100 microseconds following 
initial impact. 

« 

The significant dip in stress which occurs near 
1000 microseconds occurs when there is a change in axial 
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Figure III. 10. Nearest element stress history for near 

field analysis. / 




Figure III. 11. Nearest element effective plastic strain 

time history for near field analysis. 
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direction of the two end plates. As the relative displacement 
between the two end plates becomes zero, the stress is 
relieved. As the relative displacement continues, the stress 
once again builds. 

Similar plots are provided for the shell elements 
located 90 and 180 degrees along the circumference from the 
nearest shell element. The element at 90 degrees will be 
called the mid element and the one at 180 degrees will be 
called the remote element. 

The plots for the mid element are Figures III. 12 
and III. 13. Again, once the stress build up starts, it takes 
approximately 30 microseconds for the stress to reach yield 
stress. The mid element spends very little time at yield 
stress. As a result, the effective strain remains constant 
throughout the time of analysis. It should be noted that the 
damage mechanism is different for this element than the near 
element. The primary damage mechanism for the near element 
was the force of the shock wave impinging on the structure. 
However, as stated earlier, the water shock wave effect is 
significantly reduced by the time it reaches the mid element. 
Most of the stress at the mid element is caused by the global 
deformation of the cylinder. The cylinder tends to bend at 
the center much as a beam is loaded in one direction at the 
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Figure III. 12. 



Mid element Von Mises stress time history 
for near field analysis. 




Figure III. 13. 



Mid element effective plastic strain time 
history for .near field analysis . 
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center and in the other direction at the ends. In addition, 
the surface at the mid element is force to bulge outward by 
the movement of the near element upward toward the remote 
element. As a result, the mid element ends up in tension in 
both the axial and hoop directions. However, since the mid 
element is half way between the near and remote elements in a 
plane perpendicular to the direction of shock wave travel, the 
overall affect is minimal resulting in relatively small 
effective plastic strain. 

The final set of stress and strain plots for the 
near field analysis are for the remote element. They are 
shown in Figures III. 14 and III. 15. As in the mid element, 
the shock wave in the water at this location is mostly 
dissipated. As a result, this is not a major damage mechanism 
at the remote element. However, two other factors appear to 
be important. First, the shock wave travels circumferentially 
through the shell direction in the two semi-circular paths 
from the nearest element to the remote element because of the 
symmetry. This causes an intensification of stress at the 
most remote element. This intensified stress corresponds to 
the sharp peak in Figure III. 14. This intensification can 
also be seen in the fringe plot of Figure III. 16. 
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Figure III. 14. Remote element Von Mises stress time 

history plot for near field analysis. 




Figure III. 15. Remote element effective plastic strain 

time history plot for near field analysis. 
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The second, and probably most important damage 
mechanism for the remote element is the beam bending type 
stress placed on the cylinder. This causes a high tensile 
stress in the axial direction in the remote element while the 
flattening of the cylinder caused by the deformation at the 
center on the near side places a relatively high compressive 
sti'ess in the hoop direction. 



time - . 23200E+03 

fringes of ef f . stress ( v-m ) 

min* 2.059E-06 in element 271 

max* 4.100E-03 in element 33 



fringe 1 eve 1 s 

I 6.310E-04 
1 .341E-03 
2 . 051E-03 
2 . 762E-03 
3.472E-03 




Figure III. 16. Stress intensification in remote element 

during near field attack. 

The entire deformation occurs very rapidly and is 
so extreme that very little oscillatory response is seen in 
the structure. The only oscillatory response noted was at a 
very low frequency (approximately 500 Hz) in the axial 
direction. No whipping or breathing modes were noted. This 
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will be contrasted to the far field model later in this 
document. The final deformed shape of the cylinder is shown 
in Figure III. 17. 
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Figure III. 17. Final deformed condition of cylinder for 

near field attack. 



It should be noted that no material failure model was used 
for this analysis. The maximum effective strain was in the 
neighborhood of 18 percent. Such large strains could result 
in material failure if a failure model was considered. 
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b . STIFFENED MODEL 



Figure III. 18 shows the final deformed shape for 
a stiffened cylinder. Damage modes for the shell are very 
similar to the modes noted for the unstiffened cylinder. 
However, the stiffened cylinder withstood a larger charge at 
the same standoff distance. 

The only additional information found in this 
portion of the analysis was the failure mode of the stiffener. 
It can be clearly seen in Figure III. 18 that the stiffener 
fails by local buckling (tripping) . This is the expected mode 
of failure. This observation clearly demonstrates that it is 
important to provide lateral stiffeners to stabilize the 
circumferential stiffeners if maximum effect is to be obtained 
from the circumferential stiffeners. 
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Figure III. 18. 



Stiffener failure by tripping in near 
field attack stiffened model. 
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IV. FAR FIELD SIDE-ON EXPLOSION ANALYSIS 



A. GENERAL DESCRIPTION AND OBJECTIVES 

This study compared predictions obtained from a numerical 
analysis to the results obtained from an underwater explosion 
test. There were two main objectives. The first was to 
determine if numerical methods can adequately model the non- 
linear response of a simple cylinder subject to a side-on 
underwater explosion. If a close correlation is shown, this 
study could be used as a stepping stone to the study of more 
complex structures and materials with the final objective of 
using this type of modelling as a research tool for 
understanding the response of ships and submarines to 
underwater explosions and as a design tool for optimizing ship 
and submarine hull structures. Included in this portion of 
the study were various sensitivity analyses to determine the 
relative importance of various physical and numerical modeling 
factors on the final results. By doing this, it was hoped 
that some insight could be gained in improving future 
modelling efforts. 

The second objective was to study the non-linear response 
of the shell in an effort to obtain a better understanding of 
the possible modes of failure and response modes of the 
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cylinder. As in the near field attack, this understanding is 
necessary to be able to predict the sudden collapse of a 
cylinder subjected to underwater far field explosions. Such 
an understanding would ultimately allow the optimization of 
underwater structures to withstand this type of attack. 

B. DESCRIPTION OF MODELS USED FOR THE FAR FIELD STUDY AND 
EXPERIMENTAL METHOD 
1. PHYSICAL MODEL 

The physical model was an unstiffened right circular 
cylinder with the following characteristics. 



Dimensions : 

Length 42 inches (1.067 m) 

Diameter 12 inches (0.305 m) 

Weight 60.5 pounds (27.5 Kg) 

Materials : 



Shell 1/4 inch thick 6061-T6 Aluminum (0.64 cm) 

End Plates 1 inch thick 6061-T6 Aluminum (2.54 cm) 



The cylinders used for this test were constructed from 
commercially available material. Fabrication was performed at 
the Naval Postgraduate School. The end plates were welded to 
the shell using a Tungsten Inert Gas (TIG) process. 
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The 6061-T6 aluminum was selected on the basis of its 
high strength and strain rate insensitivity. The material 
properties of the aluminum used for the shell were verified 
using the MTS machine at the Naval Postgraduate School . 
Results of tensile testing determined that the material 
properties were near nominal with a Young's modulus of 10800 
ksi (75.6 GPa) and yield strength of 43 ksi (300 MPa). 

2. UNDERWATER EXPLOSION TEST 

The underwater explosion test was performed at the 
Dynamic Testing Incorporated (DTI) facilities in Rustburg, 
Virginia. The facility is in a quarry and the depth of the 
water is approximately 130 feet (39.6 m) at the location of 
the test. As a result, bottom reflection was not a factor in 
the test. 

The charge used for the test was 60 pounds (27.3 Kg) 
of HBX-1. The peak pressure generated by the charge was 2360 
psig (16.3 MPa) which was substantially lower than the 
calculated peak pressure of 2680 psig (18.5 MPa) for the 60 
pound (27.3 Kg) charge at a 25 foot (7.62 m) standoff 
distance. The test charge was activated by a radio control 
device . 

The test depth for both the charge and the cylinder 
was 12 feet (3.66 m) . This depth allowed the bubble generated 
by the explosion to vent at the surface prior to encountering 
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the cylinder and eliminated the possibility of a bubble pulse. 
In addition, the 12 foot (3.66 m) depth provided a clear 
pressure cutoff. 

The cylinder was held in place with a crane rig and 
the charge was suspended from a float. Distance and alignment 
of the charge to the cylinder was established and maintained 
using a tensioned span wire from the charge float to the 
cylinder support rig. Post-shot calculations found the 
arrival time of the shock wave to be consistent with a 
distance of 25 feet (3.66 m) and sound of speed in water of 
4800 ft/sec (1463 m/s). Test profile and arrangement are 
provided as Figures IV. 1 and IV. 2. 

Strain measurement was performed using CEA-06-250UW- 
350 strain gages. These are general purpose strain gages with 
an optimum range of ± 1500 microstrain and are good for both 
static and dynamic test measurements. The strain gages were 
bonded to the cylinder using a M bond 200 by a instrumentation 
technician employed by DTI. All pre-shot calibration and 
connection were performed by DTI technicians. 

The test called for 14 total strain gages (seven to 
measure hoop strains and seven to measure axial strains) . Of 
the fourteen strain gages, three failed. The dynamic range of 
the test exceeded the optimum range of the strain gages by a 
significant factor. This is the most probable cause of the 
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high strain gage failure rate. The instrumentation diagram 
for the test is provided as Figure IV. 3. The strain gage 
located at B1 was placed nearest the charge during the test. 
Strain gage output was filtered at 2000 Hz. Locations noted 
on Figure IV. 3 will be used for reference throughout the 
remainder of the report. 

Slight damage to the cylinder was noted upon 
completion of the test. Post-shot investigations found all 
strain gages firmly attached to the cylinder at the locations 
specified in the instrumentation diagram. However, some water 
intrusion was noted under the protective coating of several of 
the strain gages. This intrusion may also have played a part 
in the strain gage failures. The results of the test were 
forwarded to the Naval Postgraduate School in reference 14. 

3. NUMERICAL MODEL 

This study was performed using two primary mesh 
densities. The low density, full model mesh (Figure IV. 4) 
was used for rotation, shell type and quadrature sensitivity 
analyses. The high density quarter model was used to perform 
direct comparison to experimental results and examine end 
effects. The computational efficiency of the quarter model 
allowed a more refined mesh without a subsequent increase in 
computational time or random access memory storage capability. 
A sample quarter model was run and results checked against a 
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full model with the same mesh configuration to certify that 
the symmetry boundary conditions used to form the quarter 
model were valid. The refined mesh quarter model is shown in 
Figure IV. 5. 




Figure IV. 4. Low density, full model 

In addition to the two models noted above, several 
additional quarter models with varying mesh density were run 
to verify mesh size independence of the quarter model results. 
It was found that the most critical locations for the mesh 
sensitivity check were the locations with the highest strain. 
The areas with the highest strain were located near each end 
on the side of the cylinder located nearest the explosive 
charge. Figure IV. 6 shows the strain pattern on the surface 
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Figure IV. 5. Refined mesh model. 



on the surface of cylinder side nearest the charge. The high 
strain locations are symmetrically located 16.5 inches (0.42 
m) from the axial midpoint of the cylinder. The other region 
of significant plastic strain was located on the surface of 
the reverse side of the cylinder at the axial midpoint. 
Figure IV. 7 shows the effective plastic strain pattern for 
this location. The near side high strain regions cover a much 
smaller area than the reverse side region. That is, much 
higher strain gradients occurred on the near side compared to 
other locations on the cylinder. This condition plays a 
significant roll in mesh design and integration time increment 
selection . 
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Figure IV. 6. Effective plastic strain pattern on cylinder 
side nearest the explosive charge. 
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Figure IV. 7. Effective plastic strain pattern on cylinder 
side most remote from the charge. 
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Figures IV. 8 through IV. 10 show the results of the 
mesh sensitivity test. It was found that strains in the axial 
direction were more sensitive to mesh density than hoop strain 
results. Figure IV. 8 shows the strain at the surface of the 
cylinder at the point nearest the charge (location Bl) . This 
location has no permanent plastic strain. It can be seen that 
there is almost no significant difference between the results 
for the three mesh densities checked. Figure IV. 9 shows 
strain results for the surface of the shell at the point most 
remote from the charge in the circumferential direction at the 
axial midplane (location B3 ) . This location had the second 
highest strain of the positions checked. Although there is a 
slight difference between the three different mesh results, it 
is apparent that these differences are insignificant when 
compared to the overall plastic strain. Figure IV. 10 shows 
the strain results for the locations that experienced the 
highest plastic strains (locations A1 and Cl) . The difference 
in the hoop direction is noticeable but small enough to be 
neglected. However, the results in the axial direction are 
significant with a 30 percent variance between the average 
plastic strains for the high density mesh and medium density 
mesh. Additional refinement was not possible due to random 
access memory limitations on the system used to perform the 
analysis. On the basis of the above results it was determined 
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Figure IV. 8. Mesh sensitivity comparison for surface of 
shell located nearest the explosive charge 
(location Bl) . 
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Figure IV. 9. Mesh sensitivity comparison for position 
with largest plastic strain (locations Al 
and Cl) . 
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Figure IV. 10. Mesh sensitivity comparison for point on 

cylinder most circumerentially remote from 
the charge (location B3) . 
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that the medium mesh model was adequate for comparison of 
numerical to experimental results for all hoop strains and all 
axial strains except at the locations near the end on the side 
nearest the charge. The high density mesh was used for the 
axial strain comparison at the remaining locations. Care was 
taken to ensure that the mesh was as uniform as possible for 
both the full and the quarter mesh models to avoid problems 
with mesh reflection as noted earlier in this report. 

Thin shell elements were used for both the shell and 
end plates. Since relatively small out of plane displacements 
were encountered in the test model, it was determined that the 
four node Belytschko/Lin/Tsay shell formulation, which is the 
default formulation for VEC/DYNA3D, was adequate for the 
analysis. A Hughes/Liu (Ref. 15) shell model and a eight node 
brick shell model were also run for comparison. 

The Belytschko/Lin/Tsay shell was selected over the 
Hughes/Liu shell and 8 node brick shell formulation because of 
its higher relative computational efficiency. 

The aluminum was treated as a kinematic/isotropic 
elastic/plastic material with no strain rate sensitivity. 
Research has shown that shock velocities much higher than the 
velocities encountered in the test are required to induce 
strain rate sensitivity in 6061-T6 aluminum. 
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The pressure input for the model was obtained from the 
free field pressure transducer time record of the underwater 
explosion test. The 17000 point trace was numerically 
condensed to 100 points and entered into the TIMINT pre- 
processor of USA using the VARLIN (variable linear) option. 
Figure IV. 11 shows the pressure profile used for the analysis. 
Free surface effects were neglected and the speed of sound in 
water used for the test was 4800 ft/sec (1463 m/s) since the 
test was performed in fresh water at approximately 40 degrees 
Fahrenheit (4.5 degrees centigrade) . 




TIME (sec) 



Figure IV. 11. Undex pressure profile. 



75 



C. FAR FIELD STUDY RESULTS 

1. EXPERIMENTAL TO NUMERICAL COMPARISON 

As described earlier in the report, an underwater 
explosion test was conducted at the Dynamic Testing, 
Incorporated facility in Rustburg, Virginia. The test 
included a side-on attack of a cylinder with a stand off 
distance of 25 feet (7.62 m) using a 60 pound (27.3 Kg) HBX-1 
charge. Fourteen strain gages were attached to the cylinder, 
of which eleven provided useable data. Four statements can be 
made about the results. First, the numerical results compared 
well with the experimental results qualitatively. That is, the 
numerical response had the same general shape as the 
experimental results and it predicted compression and tension 
correctly. There was one exception to the above statement at 
position B3 (Figure IV. 20). The numerical model indicated a 
tensile axial strain at position B3 while the experimental 
data indicated a compressive strain. Physically, it can be 
observed that the shock wave is spherical and initially 
strikes the cylinder center. This places the cylinder in 
bending. Therefore, tensile strain is expected in the axial 
direction on the reverse side of the cylinder. It is believed 
that the poles on the axial strain gage at position B3 were 
reversed resulting in an error in sign of the data returned by 
the strain gage. As a result, the negative of the 
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experimental strain is plotted versus the numerical results in 
Figure IV. 20 with satisfactory results. 

Second, there were variations in magnitude between the 
numerical results and the experimental data. Further, 
magnitudes matched the experimental results more closely at 
the position nearest the charge and error increased as 
distance from the point nearest the charge increased in both 
the axial and circumferential directions. In addition, 
numerical and experimental results match more closely in areas 
with lower values of total strain. Finally, axial strains 
were affected more than hoop strains. Charge size factors 
were eliminated as a possible cause of the magnitude 
differential since the measured pressure profile was used to 
perform the post underwater explosion test numerical 
calculations. In addition, the possibility of the charge being 
located closer to the shell than the specified standoff 
distance was eliminated by comparing the actual shock wave 
travel time measured from the strain gage traces to the 
expected shock wave travel time calculated for the speed of 
sound in water for fresh water at 40 degrees Fahrenheit (4.4 
degree centigrade) . The results indicated less than two 
inches difference between the calculated and measured values 
for stand off distance. 
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Third, the frequency of oscillation of the numerical 
data was lower than the experimental data. The higher 
frequency oscillation in the physical model compared to the 
numerical model indicates that the experimental model is 
" stiff er" than the numerical model. This is an unexpected 
result, since numerical finite element solutions are normally 
expected to be stiffer than the physical model. In addition, 
the numerical results for axial strain tended to "ring" at all 
locations. The "ringing" is not a significant factor for hoop 
strains. It should also be noted that the "ringing" is 
heaviest at the front and back of the cylinder at the center. 
The causes of the "ringing" and the high relative stiffness of 
the physical model have not been determined and are a topic of 
additional study. 

Finally, there is an unexpected asymmetry in the 
experimental results. The axial strain gage at position Cl 
(Figure IV. 21) measured 50% lower than the axial strain gage 
at Al (Figure IV. 12) and the hoop strain gage at position C2 
(Figure IV. 22) measured nearly 50 percent higher than the hoop 
strain gage at position A2 (Figure IV. 13) . Failure of strain 
gages at positions A1,C1, and C2 prevented additional 
comparisons. The asymmetric results can result from two 
factors. The shell may have been rotated from the expected 
orientation by underwater currents or by forces placed on the 
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cylinder and rigging by the instrumentation cables or there 
could have been a failure in the bonding between the strain 
gage and the cylinder surface on one or more strain gages. 

Figures IV. 12 through IV. 22 provide the results of the 
numerical to experimental data comparison. 
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Figure IV. 12. Experimental /numerical comparison for 

position Al axial strain. 
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IV. 13. Experimental /numerical comparison for 
position A2 hoop strain. 
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Figure IV. 14. Experimental /numerical comparison for 

position A2 axial strain. 
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Figure IV. 15. Experimental/numerical comparison for 

position B1 hoop strain. 
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Figure IV. 16. Experimental/numerical comparison for 

position B1 axial strain. 
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IV. 17. Experimental/numerical comparison for 
position B2 hoop strain. 
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Figure IV. 18. Experimental/numerical comparison for 

position B2 axial strain. 
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Figure IV. 19. Experimental/numerical comparison for 

position B3 hoop strain. 
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Figure IV. 20. Experimental /numerical comparison for 

position B3 axial strain. 
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Figure IV. 21. Experimental /numerical comparison for 

position Cl axial strain. 
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Figure IV. 22. Experimental /numerical comparison for 

position C2 hoop strain. 
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2. SENSITIVITY ANALYSES 



A series of sensitivity analyses were performed in an 
effort to explain the differences between the numerical and 
experimental results noted in the previous section. In 
addition, these analyses provided additional insight into the 
relative importance of various factors in the performance of 
underwater explosion tests and the associated calculations. 
Seven sensitivity analyses were performed. The first was the 
mesh sensitivity test. The results of this analysis have 
already been discussed. The other six analyses were, end 
effect, shell element formulation, integration time increment 
length, quadrature, rotational position and r| value 
sensitivity checks. The results of these analyses are 
provided in the following subsections. 

a , END EFFECT SENSITIVITY ANALYSIS 

As previously noted, the most severe deformation 
occurred at locations near the end of the cylinder (positions 
A1 and Cl) . Two processes cause this phenomena. First, the 
relatively large mass of the end plates apply large inertial 
forces to the cylinder shell near the end plates. Second, the 
one inch thick end plates are very stiff and their lack of 
flexibility causes the weaker material of the shell near the 
end plates to deform in response to applied forces. A 
examination of the numerical and experimental data reveals 
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that these effects are concentrated near the end plates and 
result in large strain gradients. This means that elements on 
either side of a selected element near the end of the cylinder 
can have significantly different strain values. Accurate 
placement of strain gages within this region and careful mesh 
design along with adequately short time integration increments 
are critical in obtaining satisfactory results in a numerical 
to experimental data comparison. In addition, as stated 
earlier, the end plates are attached to the shell using a 
Tungsten Inert Gas process. This welding process results in 
high temperatures near the end of the cylinder. Since the 
aluminum for this model is at a peak hardened condition, this 
process could result in a change of the material properties 
near the end of the cylinder that can only be restored by 
performing the age hardening process again after the welding 
is complete. These factors can result in an uncertainty in 
the expected strain compared to what might occur under ideal 
circumstances . 

The mesh sensitivity results clearly display the 
importance of mesh design within this region. However, even 
with proper design, the large gradients can result in 
significant differences between the predicted and actual 
strains since the strain computed for the element is an 
average of the strain over the entire element vice a strain at 
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a specific point. The best possible results would be obtained 
in these regions with large gradients if the mesh could be 
refined such that the size of the elements is the same size as 
the gage length of the strain gage. However, this would 
result in a prohibitively large number of elements and a 
subsequent increase in problem solution times. These problems 
can be overcome by placing strain gages in areas that are 
expected to have consistently increasing or decreasing strains 
and then ensuring that the mesh is designed so that the strain 
gage location is at the center of the element. If possible, 
large gradient regions should be avoided. If strain gages 
must be placed in a high gradient region, then the strain 
gages should be placed to one side or the other of the minimum 
or maximum strain location. Placement at the minimum or 
maximum point will result in an error since the average for 
the element will lie above a minimum or below a maximum if the 
element is not the same size as the gage length of the strain 
gage . 

In this study, the strain gages located at Al and 
Cl were located at the point of highest compressive strain. 
Therefore, a study was performed to determine the relative 
importance of the noted location factors. Figures IV. 23 
through IV. 27 show the results of this study. Strains of two 
additional elements nearer the end were compared to the 
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Figure IV. 24. End effect sensitivity results. 
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Figure IV. 25. End effect sensitivity results. 
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End effect sensitivity results. 
(Cl Axial) 



89 




TIME (SEC) 

Figure IV. 27. End effect sensitivity results. 
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measured strain and the actual strain gage location. Elements 
nearer the end plate were selected since the welding effects 
described in the previous paragraph would tend to move the 
high strain location nearer to the end plate by weakening the 
material near the end plate. Only the positions with useable 
experimental results are shown. In four of the five cases 
(positions A1 axial, A2 hoop, Cl axial and C2 hoop), if 
asymmetry effects are taken into account, the element one 
nearer to the end from the actual strain location provides a 
better estimate of' the actual strain measured during the 
underwater explosion test. At the fifth location (A2 axial), 
the second element closer to the end provides the best 
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results. These results require additional study to separate 
and quantify the effect of the phenomena. 

b. SHELL FORMULATION, QUADRATURE RULE AND INTEGRATION 

TIME INCREMENT SENSITIVITY ANALYSES 

In addition to the above end effects, there was 
some concern that the mid plane reference for the thin shell 
element would result in a greater flexible length than the 
actual physical model . This concern was based on the fact 
that the mass and stiffness of the end plates is concentrated 
into a planar surface co-located with the mid plane of the end 
plate in the thin shell analysis. This resulted in the shell 
portion of the structure being one inch longer in the 
numerical model than the physical model . This problem could 
have been avoided by using the Hughes/Liu formulation and 
shifting the reference plane to the inner surface of the 
shell. To resolve this issue a study was conducted to compare 
the performance of different types of thin and thick shell 
element formulations. 

Results from the Belytschko/Tsay/Lin shell 
formulation were compared to results from the same model using 

the Hughes/Liu shell formulation. As stated earlier, the 

< 

Belytschko/Lin/Tsay shell has the advantages of increased 
computational efficiency and a high degree of stability with 
large deformations at the expense of reduced accuracy at high 
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levels of plastic strain. The major difference between the 
two formulation stems from the fact that the element normal 
direction is updated periodically in the Hughes/Liu 
formulation. The Belytschko/Lin/Tsay formulation assumes 
negligible out of plane deformations, and therefore, does not 
update the shell normal. As a result, the inaccuracy of the 
Belytschko/Lin/Tsay formulation will increase as shell 
deformation becomes significant. 

The models used to compare the two formulations 
were identical in all aspects with the exception of the shell 
formulation. The center line plane was used for the reference 
on both models. The results confirmed that the strain levels 
encountered in this underwater explosion test were small 
enough to support use of the Belytschko/Lin/Tsay formulation. 
However, it was apparent that differences did occur for 
positions with significant plastic strain in the axial 
direction (Positions Al, A2 , B3, Cl, and C2 ) . Although the 
differences in these cases were not significant enough to 
require use of the Hughes/Liu formulation, it is also noted 
that higher strain may result in larger differences. 
Therefore the Belytschko/Lin/Tsay formulation should not be 
used in cases where significant denting occurs unless 
stability problems occur while using the Hughes/Liu 
formulation . 
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As stated earlier, the presence of high strain 
gradients near the end plates causes small changes in end 
condition or distance to be significant. When it became 
apparent that end effects would be important in the results an 
investigation was performed to determine if an eight node 
brick shell formulation would provide more accurate results 
near the end of the cylinder. The thin shell formulation 
results as well as the experimental results were compared to 
results from a model computed using eight node brick shell 
elements. All three formulations are compared to experimental 
results in Figures IV. 28 through IV. 38. The following 
information can be gleaned from the plotted results. First, 
it is apparent that the greatest differences occur near the 
positions with the highest strains. At the same time, it can 
be noted that there is virtually no difference at the 
locations with no permanent strain. Second, as shown in 
Figures IV. 28, IV. 29, IV. 35, IV. 36 and IV. 38, it is clear that 
there is a significant difference between the eight node brick 
shell results and the Belytschko/Lin/Tsay results at the 
locations with high levels of permanent plastic strain. 
However, contrary to the expected results, the eight node 
brick shell results move further from the expected values than 
the other formulations. It is also noted that the Hughes/Liu 
formulation lies between the eight node brick shell and the 
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Figure IV. 28. Shell formulation sensitivity results. 
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Figure IV. 29. Shell formulation sensitivity results. 
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Figure IV. 30. Shell formulation sensitivity results 
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Figure IV. 31. Shell formulation sensitivity results. 
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Figure IV. 32. Shell formulation sensitivity results. 
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Figure IV. 33. Shell formulation sensitivity results. 
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Figure IV. 34. Shell formulation sensitivity results. 
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Figure IV. 35. Shell formulation sensitivity results. 
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Figure IV. 36. Shell formulation sensitivity results. 
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Figure IV. 37. Shell formulation sensitivity results. 
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Figure IV. 38. Shell formulation sensitivity results. 
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Belytschko./Lin/Tsay formulation . 

Additional research was performed to determine the 
cause of the disparities. The study revealed that the eight 
node brick shell is sensitive to integration time increment 
and will move marginally closer to the thin shell results if 
time integration is cut in half. However, the overall shift 
is only about 10 percent of the total difference. Quadrature 
rule (numbei of points used in the Gauss quadrature numerical 
integration scheme) proved to be a more significant effect. 
Use of five point quadrature moved the three results closer 
together while having the most profound effect on the 
Belytschko/Lin/Tsay formulation. Again, the affect of 
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quadrature rule affected the thick shell results only 
marginally. Figure IV. 39 shows the combined results for the 
location with greatest plastic deformation. 




Figure IV. 39. Effect of changing quadrature rule and time 

integration increment at location of highest 
strain (Al and Cl hoop) . 



In summary, all three formulations appear to be 
satisfactory as long as care is used in designing the mesh and 
selecting the integration time and quadrature integration 
rule parameters . Specifically, when using Belytschko/Lin/Tsay 
formulations in areas with relatively high strain, the number 
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of quadrature points should be increased until stable results 
are achieved. When using eight node brick shell elements, 
integration time increment must be selected with care but 
number of quadrature points seems to be less critical . The 
Hughes/Liu formulation appeared to be relatively insensitive 
to both quadrature rule and integration time increment. 

Reference 16 provides some useful thumbrules for 
selection of time increments. The following criteria are 
recommended . 

At « 0.9-r— — r- for brick shells 
(A b c) 

A 

At * 0.9— - for thin shells 

Dc 

V - element volume 
At - time increment 
A e - maximum surface area 
D - maximum diagonal 
c - speed of sound in the material 
A b - maximum area of any surface 



The above criteria were found to be adequate except for the 
highest strain areas where the thick shell element rule did 
not provide stable solutions. In areas such as Al and Cl, a 
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value of the integration time increment half of the above 
recommendation proved to be satisfactory for the eight node 
brick shell. 

C. ROTATION SENSITIVITY ANALYSIS 

A sensitivity analysis was performed to determine 
the effect of an in plane rotation away from the expected 
symmetric orientation in an effort to explain the cause of the 
asymmetric results of the underwater explosion test. It was 
hypothesized that an unplanned rotation greater than ten 
degrees would have been detected by the personnel performing 
the test. Four different models were run within this range 
representing rotations of 0.0, 2.5, 5.0 and 10.0 degrees. The 
results are shown along with experimental results where 
available in Figures IV. 40 through IV. 53. The following 
observations are made concerning the results. First, the most 
dramatic affects are on the reverse side of the cylinder at 
position B3 (Figures IV. 48 and IV. 49) . The results show that 
the differential between the numerical and experimental 
results at position B3 can be explained by a six to eight 
degree rotation from the symmetric configuration. Rotational 
effects at locations B1 and B2 on the centerline (Figures 
IV. 44 through IV. 47) are insignificant. Hoop strain at 
position C2 (Figure IV. 52) is approximately 60 percent 
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Figure IV. 40. Rotation sensitivity results. 
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Figure IV. 41. Rotation sensitivity results. 
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Figure IV. 42. Rotation sensitivity results. 
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Figure IV. 43. Rotation sensitivity results. 
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Figure IV. 44. Rotation sensitivity results. 
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Figure IV. 45. Rotation sensitivity results. 
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Figure IV. 46. Rotation sensitivity results. 
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Figure IV. 47. Rotation sensitivity results. 
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Figure IV. 48. Rotation sensitivity results. 
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Figure IV. 49. Rotation sensitivity results. 

(B3 Axial) 
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Figure IV. 50. Rotation sensitivity results. 

(Cl Hoop) 




TIME (SEC) 



Figure IV. 51. Rotation sensitivity results. 
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Figure IV. 52. Rotation sensitivity results. 
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Figure IV. 53. Rotation sensitivity results. 

(C2 Axial) 
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higher than the hoop strain at position A2 (Figure IV. 43) 
with a rotation of ten degrees. This is also consistent with 
the experimental data. Similar positive results were obtained 
for positions A1 and Cl axial strains. It was further 
discovered that rotating the cylinder about its axis could 
further improve the results. However, even though these 
rotations did improve the results, significant differences 
still exist between the experimental and numerical strains at 
the ends of the cylinders. Although it is clear that the 
model can account for rotational effects, it is also clear 
that other factors are causing the large differences. Once 
again, welding affects are suspected to be the probable cause. 

The important point to note out of these results is 
that even small rotations from expected orientation can result 
in significant errors on in expected results. Therefore 
extreme care must be taken to ensure that instrumentation 
cable tension and other unanticipated factors do not cause 
undetected rotations. 

3 . PHYSICAL FINDINGS 
a . RESPONSE MODES 

It was determined that a cylinder subject to a side 
on explosion will have three primary response modes. The 
first mode is an accordion motion. The accordion motion 
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results from the compression and subsequent release of the 
cylinder in the axial direction. Figure IV. 54 shows a plot of 
points located at the center of each end plate. It is clear 




TIME (sec) 

Figure IV. 54. Cylinder accordion motion. 

that the two end plates are travelling in opposite directions 
at the same time generating the accordion motion. 

The cylinder is also subject to a whipping mode 
parallel to the direction of shock wave travel. The whipping 
mode is the most significant motion experienced by the 
cylinder and is caused as a result of the curvature of the 
shock wave. In the symmetric situation, the shock wave will 
come in contact with the center of the cylinder first. This 
will cause the center to move first, followed by the ends. 
The cylinder will then move in an oscillatory motion that is 
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a function of the stiffness and mass distribution of the 
cylinder and the water surrounding the cylinder. Figure IV. 55 




TIME ( sec ) 

Figure IV. 55. Cylinder whipping motion in plane parallel 

to shock wave direction. 

shows a plot of a points located at the center and ends of a 
line located parallel to the axis on the near side of the 
cylinder. The plot shows that the end plates are moving in 
the opposite direction of the cylinder throughout the 
transient response of the cylinder. Figure IV. 56 shows a 
scale factor 20 drawing of the cylinder at two different 
times. The cylinder's opposite direction of curvature at the 
two different times is a result of the whipping motion. 

The final response mode noted was a breathing 
motion in the plane perpendicular to the shock wave direction 
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Figure IV. 56. Cylinder curvature as a result of whipping 

motion (scale factor 20). 

of travel. Although breathing motion also occurred in the 
direction parallel to the shock wave travel, it was not as 
obvious since the much larger whipping motion turned out to be 
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the predominant mode in that direction. Figure IV. 57 shows a 
plot of two points located at the top and bottom of the 




Figure IV. 57. Cylinder breathing motion perpendicular to 

the shock wave direction of travel. 

cylinder in a plane perpendicular to the axis at the axial mid 
point of the cylinder. It can be observed that the upper 
point is moving in a direction opposite to the lower point 
throughout the transient response of the cylinder. The 
breathing motion is also caused by the compression and 
subsequent release of the cylinder. Figure IV. 58 provides an 
illustration of the breathing motion. The two scale factor 40 
drawings are for two separate times and show the shell first 
bowed inward toward the axis and then outward away from the 
axis . 
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Figure IV. 58. Illustration of cylinder breathing mode at 

two different times (scale factor 40). 

Jt>. ROTATIONAL EFFECTS 

Plastic strain fringe plots generated as a result 
of the rotation sensitivity analysis revealed some interesting 
information on the causes of the strain distribution 
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experienced by the cylinder. The experimental results 
included a reduction in the strain at the rear of the cylinder 
at position B3 , a decrease at A2 relative to C2 and an 
increase at Al relative to Cl. The fringe plots show why this 
strain distribution occurs. Figure IV. 59 shows the effective 
plastic strain distribution for a 7.5 degree rotation. The 
left side of the cylinder is nearest the charge. The results 
show that the rotation tends to diffuse the strain around the 
cylinder on the near end while concentrating it at the far 
end. This causes the distribution noted for positions Al , Cl, 
A2 and C2 . At the same time, the high stress region on the 
reverse side of the cylinder tends to move away from the 
charge. This placed location B3 in a lower strain region 
which led to the experimental and numerical results noted at 
position B3 . 
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Figure IV. 59. Effective plastic strain distribution for 

near and remote side of cylinder with 7.5 
degrees of rotation. 
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V. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS. 

1. NUMERICAL MODELLING 

Two general conclusions can be reached from the 
material contained in this report. First, the USA/DYNA3d 
connection is successful and can replicate the response of 
simple analytical models. 

Second, numerical modeling can predict the response of 
a simple cylinder to an underwater explosion. Near field 
numerical predictions matched expected results. However, to 
fully confirm this, experimental data must be obtained for 
comparison to the numerical results. 

Far field numerical predictions generally match 
experimental results if rotation and end effects resulting 
from fabrication caused material property changes are 
correctly modelled. It was found that results in high strain 
areas are extremely sensitive to shell formulation, mesh 
design, quadrature rule and integration time increment. The 
best results were achieve with brick shell elements. However, 
the eight node brick shell required substantially longer 
computation times to achieve the desired results because of 
the need to reduce the integration time increment. In 
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addition, it was found that thin shell formulations can also 
provide correct results. However, results for the 
Belytschko/Lin/Tsay formulation appear to be very sensitive to 
the number of quadrature points used for the numerical 
integration scheme in high strain areas. 

2. PHYSICAL ASPECTS 

a. NEAR FIELD PHYSICAL RESULTS 

The propagation of the pressure wave through the 
water near the shell was studied. The numerical model 
correctly predicted pressure response at the shell water 
interface. In addition, cavitation was predicted next to the 
shell . 

Damage to the shell was primarily local in nature 
with severe buckling occurring at the point nearest the 
charge. In addition, stiffeners in stiffened models failed by 
local buckling (tripping) . Further the remote side of the 
shell experience significant plastic strain as a result of the 
bending effect caused by the shock wave striking the axial mid 
plane of the cylinder. 

b. FAR FIELD PHYSICAL RESULTS 

Primary damage areas are near the ends of the 
cylinder on the side nearest the charge where the stiff, 
heavy, flat end plates caused a concentration of the effective 
plastic strain. Damage also occurred on the reverse side as 
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a result of a bending effect similar to that described in the 
near field results. The cylinder experienced breathing, 
whipping and accordion response modes. 

In addition, it was discovered that rotation tends 
to diffuse strain on the end nearest the charge while 
concentrating the strain at the far end on the side nearest 
the charge. The high strain area located at the center of the 
cylinder on the reverse side tends to migrate toward the end 
most remote from the charge. 

B . RECOMMENDATIONS . 

1. TOPICS FOR ADDITIONAL STUDY. 

a. WELDING FABRICATION EFFECTS. 

An analysis should be performed to quantify the 
relative effect that the change in material properties 



genei ated 


by the 


welding fabrication 


process 


has on 


the 


numerical 


results 


This analysis 


could 


include 


the 


measurement of material properties near 


a weldment . These 



properties could then be added as a separate material in the 
numerical model. 

b. EIGHT NODE BRICK SHELL SENSITIVITY ANALYSES. 

Although it was fairly clear that the eight node 
brick shell formulation comes closest to predicting the 
overall response of the shell, it was also noted that the 
formulation is very sensitive to integration time increment in 
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areas with high strain. Commonly used thumbrules did not 
appear to be adequate in this case. In addition, additional 
analyses need to be performed to determine the mesh 
sensitivity of the eight node brick shell in this model . 

C. FAILURE CRITERIA. 

This study was performed on a model with 
relatively low total plastic strain (less than one percent). 
In order to deal with larger strains, a failure model must be 
introduced into the material modelling of the cylinder. The 
model should include structural instability as well as 
material rupture criteria. 

d. NEAR FIELD EXPERIMENTATION 

Although the numerical predictions appear to be 
physically correct, the physical results obtained using them 
cannot be assumed completely correct until they are confirmed 
with experimental results. A study should be conducted to 
compare near field experimental results with numerical 
predictions . 

2. RECOMMENDATIONS TO IMPROVE TEST CONTROL. 

Several factors made the comparison of the numerical 
to experimental results difficult. If properly controlled, 
the analysis process could be simplified. First, rotation of 
the cylinder must be carefully controlled. Second, unless 
specifically required, high strain gradient areas should be 
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avoided. Placement of the strain gages becomes critical in 
these locations as does mesh design and integration time 
increment. If these areas cannot be avoided, additional 
sensitivity analyses may be required to determine the adequacy 
of the mesh and integration time increment. Finally, analysis 
near welded seams should be avoided unless the effects can be 
quantified. If near weld analysis cannot be avoided, 
consideration should be given to restoring the heat treatment 
after the weld process. 
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